{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "reinforce.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyPsMNqYXtbNkyFhna612LEV",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/letianzj/QuantResearch/blob/master/ml/reinforce.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JlNkxcEuDTVn",
        "colab_type": "text"
      },
      "source": [
        "## Monte Carlo Policy Gradient CartPole\n",
        "\n",
        "Use MC-PG, a.k.a REINFORCE, to solve CartPole game.\n",
        "\n",
        "[OpenAI Gym CartPole](https://gym.openai.com/envs/CartPole-v0/) has four states, cart position, cart speed, pole angle, and pole speed. The actions are either going left or right. The objective is to keep pole from falling. Every move that doesn't lead to a fall gets reward 1.\n",
        "\n",
        "Monte-Carlo policy gradient is on-policy model-free sampling method. The expected total rewards,\n",
        "\n",
        "$$\n",
        "J(\\theta)=E_{\\tau \\sim \\pi_{\\theta}(\\tau)}[R(\\tau)]=\\int R(\\tau) p_{\\theta}(\\tau)d\\tau\n",
        "$$\n",
        "\n",
        "and its derivative on policy parameter $\\theta$ is,\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "\\frac{\\partial J}{\\partial \\theta}&=\\int \\left( \\frac{\\partial}{\\partial \\theta} p_{\\theta}(\\tau) \\right) R(\\tau)d\\tau \\\\\n",
        "&=\\int p_{\\theta}(\\tau) \\left(\\frac{\\partial}{\\partial\\theta}log p_{\\theta}(\\tau)\\right) R(\\tau)d\\tau \\\\\n",
        "&=E_{\\tau \\sim \\pi_{\\theta}(\\tau)} \\left[\\frac{\\partial}{\\partial\\theta}log p_{\\theta}(\\tau) R(\\tau) \\right] \\\\\n",
        "&=E_{\\tau \\sim \\pi_{\\theta}(\\tau)} \\left[\\sum_{t=1}^{T} \\frac{\\partial}{\\partial\\theta}log \\left( \\pi_{\\theta}(a_t|s_t) \\right)R(\\tau) \\right] \\\\\n",
        "&\\approx \\frac{1}{N} \\sum_{n=1}^N \\left( \\left( \\sum_{t=1}^{T} \\frac{\\partial}{\\partial \\theta} log \\pi_{\\theta} \\left( a_t^{(n)} | s_t^{(n)}\\right) \\right) R\\left(\\tau^{(n)}\\right) \\right)\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "Policy gradient updates parameters iteratively along every sample trajectory. We can calculate the gradient directly, or construct a cross-entropy loss function whose gradient is equal to the policy gradient as well explained [here](https://adventuresinmachinelearning.com/policy-gradient-tensorflow-2/). \n",
        "\n",
        "$$\n",
        "Loss= CE = \\left( \\sum_{t=1}^{T} logp_{\\theta} \\left( a_t^{(n)} | s_t^{(n)}\\right) \\right) R\\left(\\tau^{(n)}\\right)\n",
        "$$\n",
        "\n",
        "Below it uses the former and code structure follows closely [here](https://github.com/dragen1860/Deep-Learning-with-TensorFlow-book/blob/master/ch14-%E5%BC%BA%E5%8C%96%E5%AD%A6%E4%B9%A0/REINFORCE_tf.py).\n",
        "\n",
        "REINFORCE is usually adjusted by a baseline to reduce variance.\n",
        "\n",
        "\n",
        "__Reference__\n",
        "* Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.\n",
        "* [RL by David Silver](https://www.davidsilver.uk/teaching/)\n",
        "* [OpenAI Spinning Up](https://spinningup.openai.com/en/latest/)\n",
        "* [Lil'Log](https://lilianweng.github.io/lil-log/)\n",
        "* [Deep Learning with TensorFlow](https://github.com/dragen1860/Deep-Learning-with-TensorFlow-book)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0lOf7OSWHOvv",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 238
        },
        "outputId": "4ef1ee3f-33e1-48df-ea19-3eb203d72657"
      },
      "source": [
        "import gym\n",
        "import numpy as np\n",
        "from \tmatplotlib import pyplot as plt\n",
        "import tensorflow as tf\n",
        "from tensorflow import keras\n",
        "from tensorflow.keras import layers, optimizers\n",
        "\n",
        "env = gym.make('CartPole-v1')\n",
        "print(f'CartPole state {env.observation_space.shape[0]},  action {env.action_space.n}') # state size = 4, action = 0 or 1\n",
        "env.seed(5555)\n",
        "tf.random.set_seed(5555)\n",
        "np.random.seed(5555)\n",
        "\n",
        "epoch = 300\n",
        "learning_rate = 0.0002\n",
        "gamma = 0.98\n",
        "reward_to_go = True\n",
        "\n",
        "class Policy(keras.Model):\n",
        "    def __init__(self):\n",
        "        super(Policy, self).__init__()\n",
        "        self.fc1 = layers.Dense(128, kernel_initializer='he_normal', activation='tanh')\n",
        "        self.fc2 = layers.Dense(2, kernel_initializer='he_normal', activation=None)\n",
        "        self.optimizer = optimizers.Adam(lr=learning_rate)\n",
        "\n",
        "    def call(self, inputs, training=None):\n",
        "        x = self.fc1(inputs)          # [b, 4]  ==> [b 32]\n",
        "        x = self.fc2(x)               # [b, 32] ==> [b, 2]\n",
        "        # x = tf.nn.softmax(x, axis=1)       # return prob instead of logits\n",
        "        return x\n",
        "\n",
        "model = Policy()\n",
        "model(tf.random.normal((10,4)))           # batch=10, output=(10,2)\n",
        "model.summary()"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "CartPole state 4,  action 2\n",
            "Model: \"policy_7\"\n",
            "_________________________________________________________________\n",
            "Layer (type)                 Output Shape              Param #   \n",
            "=================================================================\n",
            "dense_14 (Dense)             multiple                  640       \n",
            "_________________________________________________________________\n",
            "dense_15 (Dense)             multiple                  258       \n",
            "=================================================================\n",
            "Total params: 898\n",
            "Trainable params: 898\n",
            "Non-trainable params: 0\n",
            "_________________________________________________________________\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lAmz6UWC_mex",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "scores = []\n",
        "\n",
        "for n_epi in range(epoch):             # epoch\n",
        "    with tf.GradientTape(persistent=True) as tape:\n",
        "        s = env.reset()          # return s0\n",
        "        done = False\n",
        "        log_probs = []\n",
        "        rs = []\n",
        "\n",
        "        # sample one game, one trajectory\n",
        "        while not done:\n",
        "            s = tf.constant(s,dtype=tf.float32)\n",
        "            s = tf.expand_dims(s, axis=0)                   # [4] => [1,4]\n",
        "            logits = model(s)                                    # shape = [1, 2]\n",
        "            a = tf.random.categorical(logits, 1)[0]           # categorical accepts logits, not prob\n",
        "            a = int(a)                                      # Tensor to int\n",
        "            s_prime, r, done, info = env.step(a)\n",
        "            rs.append(r)                                 # rewards\n",
        "            prob = tf.nn.softmax(logits)                 # take softmax\n",
        "            log_probs.append(tf.math.log(prob[0][a]))    # log(prob)\n",
        "            s = s_prime\n",
        "        \n",
        "        scores.append(np.sum(rs))\n",
        "        # print(f\"# of episode :{n_epi}, avg score : {scores[-1]}\")\n",
        "\n",
        "        # after game over, for each time step, update parameters\n",
        "        rs = tf.convert_to_tensor(rs, dtype=tf.float32)\n",
        "        log_probs = tf.convert_to_tensor(log_probs, dtype=tf.float32)\n",
        "\n",
        "        R = 0                                     # reward to go\n",
        "        LP = 0                                    # log\\pi to go\n",
        "        for t in range(len(rs)-1, -1, -1):\n",
        "            R = rs[t] + gamma * R\n",
        "            LP = LP + log_probs[t]\n",
        "            loss = - LP * R                 # negative sign means gradient ascent\n",
        "            with tape.stop_recording():\n",
        "                grads = tape.gradient(loss, model.trainable_variables)\n",
        "                model.optimizer.apply_gradients(zip(grads, model.trainable_variables))  \n",
        "        # rts = [np.sum(rs[i:]) for i in range(len(rs))]\n",
        "        # rts = tf.convert_to_tensor(rts, dtype=tf.float32)\n",
        "        # log_probs = tf.convert_to_tensor(log_probs, dtype=tf.float32)\n",
        "        # if reward_to_go:\n",
        "        #     loss = - tf.math.reduce_sum(log_probs * rts)         # use r(tau)_t:T; negative sign means gradient ascent\n",
        "        # else:\n",
        "        #     loss = - tf.math.reduce_sum(log_probs * rts[0])      # use r(tau); negative sign means gradient ascent\n",
        "        # grads = tape.gradient(loss, model.trainable_variables)\n",
        "        # model.optimizer.apply_gradients(zip(grads, model.trainable_variables))\n",
        "\n",
        "    del tape\n",
        "    \n",
        "env.close()"
      ],
      "execution_count": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NSz4Nd4iERFR",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "outputId": "63474ba8-6c52-4fb7-c87f-6af41da8f3bb"
      },
      "source": [
        "plt.plot(np.arange(epoch), scores)\n",
        "plt.plot(np.arange(epoch), scores, 's')\n",
        "plt.xlabel('epoch')\n",
        "plt.ylabel('score')\n",
        "plt.show()"
      ],
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5zcdX3v8ddnL9nNBTYk2UROAAM1eGmVW7QolapUK9Rj0AcqB6vUB23OsdDq8TSVam3TG9WiUm0pmhY92IpirRx4WLTFAHJ6LEhA5CJoAoaSFEgCyYbcNjs7n/PH7/ubzGzmuju/y+zv/Xw89rEzv/nNzOe3M/v7/L53c3dEREQA+rIOQERE8kNJQUREKpQURESkQklBREQqlBRERKRiIOsAZmLJkiW+YsWKrMMQEekp99577053H633WE8nhRUrVrBx48aswxAR6Slm9kSjx1R9JCIiFUoKIiJSoaQgIiIVSgoiIlKhpCAiIhWJ9j4ysy3A88AkUHL3VWa2CLgBWAFsAd7p7rvMzIDPAOcB+4Ffc/f7koxPRKa4ciXs237k9vlLYe2m9OOR1KVRUni9u5/q7qvC/cuBDe6+EtgQ7gOcC6wMP2uAa1KITUSq1UsIzbbLrJNF9dFq4Lpw+zrg/KrtX/LIXcBCMzs2g/hERAor6aTgwL+a2b1mtiZsW+buT4XbTwPLwu3lwJNVz90attUwszVmttHMNu7YsSOpuEVECinpEc2/4O7bzGwpcKuZPVr9oLu7mXW0yo+7rwfWA6xatUorBImIdFGiJQV33xZ+bwduBF4FPBNXC4XfcWXlNuD4qqcfF7aJiEhKEksKZjbfzI6KbwNvAh4CbgYuDrtdDNwUbt8MvNciZwJjVdVMIpKG+Us72y6zTpLVR8uAG6OepgwA17v7t83sHuBrZnYJ8ATwzrD/LUTdUTcTdUl9X4KxiUg9odvpSz72LVafspxPXPCKjAOStCWWFNz9ceCUOtufBc6ps92BS5OKR0TaVy5D2dVkV0Qa0SwiRyi7U1ZOKCQlBRE5QtkdV0mhkJQUROQIZVf1UVEpKYhIjXKoN1L1UTEpKYhIjbiEoJJCMSkpiEiNuISgnFBMSgoiUkMlhWJTUhCRGkoKxaakICI14uojNTQXk5KCiNSISwgap1BMSgoiUsPL0W+VFIpJSUFEakyqTaHQlBREpMbhhuaMA5FMKCmISA21KRSbkoKI1ChX2hSUFIpISUFEalSqj8oZByKZUFIQkRqV6iNUUigiJQURqeEavFZoSgoiUmOyrIbmIlNSEJEa6pJabEoKIlLj8NxHygpFpKQgIjVcJYVCU1IQkRqTGrxWaEoKIlJDg9eKTUlBRGpo8FqxKSmISA2tvFZsSgoiUiNuYFZOKCYlBRGpoZJCsQ1kHYCI5MCVK2HfdgBOB7YMA3uAK5fC2k1ZRiYpU0lBRCoJoe3tMmspKYiISIWSgoiIVCSeFMys38x+YGbfDPdPNLO7zWyzmd1gZnPC9qFwf3N4fEXSsYmISK00SgofAB6puv8J4Cp3fxGwC7gkbL8E2BW2XxX2ExGRFCWaFMzsOOBXgL8L9w14A/D1sMt1wPnh9upwn/D4OWF/EUna/KWdbZdZK+kuqX8J/C5wVLi/GNjt7qVwfyuwPNxeDjwJ4O4lMxsL+++sfkEzWwOsATjhhBMSDV6kMNZuYs/BCV6x7l9ZdvQQz+wZZ9nRQ9y99peyjkxSllhJwczeAmx393u7+bruvt7dV7n7qtHR0W6+tEihlSa95remzi6mJEsKZwFvNbPzgGHgaOAzwEIzGwilheOAbWH/bcDxwFYzGwBGgGcTjE9EqpQmoxnwJsJvTZ1dTImVFNz999z9OHdfAVwI3Obu7wZuBy4Iu10M3BRu3xzuEx6/zfWtFEnNRCgaTKikUGhZjFP4MPAhM9tM1GZwbdh+LbA4bP8QcHkGsYkU1tSSguY+KqZU5j5y9zuAO8Ltx4FX1dnnIPCONOIRkSPFJYRSOV5PQUmhiDSiWUQAKE1ZVUcFhWJSUhAR4HCvo5iqj4pJSUFEgMNtCTHVHhWTkoKIAIfbEmIqKRSTkoKIAEeWFJQTiklJQUSAw72PYiopFJOSgogAh8cpxJQUiklJQUSAeiWFjAKRTCkpiAhw5DgF0PxHRaSkICLAkeMUQKWFIlJSEBHgyN5HoHaFIlJSEBHgyHEKoKRQREoKIgIc2fsINFahiJQURAQ4svcRqKRQREoKIgLU732khubiUVIQEUAlBYkoKYgIUL9LqnJC8SgpiAigwWsSUVIQEaBR9VEGgUimlBREBKjfJVVtCsWjpCAigAavSURJQUSA+tNcKCcUj5KCiACNJsRTVigaJQURAWBCg9cEJQURCeqWFJQVCkdJQUSARuMUMghEMqWkICKAprmQiJKCiAAapyARJQURARqNU8ggEMmUkoKIAI3GKSgrFI2SgogAjcYpZBCIZEpJQUQAmNA0F0KCScHMhs3s+2b2QzN72Mz+KGw/0czuNrPNZnaDmc0J24fC/c3h8RVJxSYiR1JDs0CyJYVx4A3ufgpwKvBmMzsT+ARwlbu/CNgFXBL2vwTYFbZfFfYTkZRokR2BBJOCR/aGu4Phx4E3AF8P268Dzg+3V4f7hMfPMTNLKj4RqVWvoVklheIZSPLFzawfuBd4EXA18Biw291LYZetwPJweznwJIC7l8xsDFgM7JzymmuANQAnnHBCkuGLFIrmPsrIlSth3/Yjt89fCms3pR5Oog3N7j7p7qcCxwGvAl7Shddc7+6r3H3V6OjojGMUkUhp0hnsry2cq6SQgnoJodn2hKXS+8jddwO3A68GFppZXEI5DtgWbm8DjgcIj48Az6YRn4hE01wMD/TXbNM4heJJsvfRqJktDLfnAm8EHiFKDheE3S4Gbgq3bw73CY/f5vpGiqSmVC4zNFh7SlD1UfEk2aZwLHBdaFfoA77m7t80sx8BXzWzPwV+AFwb9r8W+Hsz2ww8B1yYYGwiMkVp0pk/p/aUoKmziyexpODuDwCn1dn+OFH7wtTtB4F3JBWPiDQQGjofMuAgMBxt3uEjbPYfZBmZZKDt6iMzm2tmL04yGBHJQIMGzVEbU5tCGuYv7Wx7wtoqKZjZfwU+CcwBTjSzU4E/dve3JhmciGRLtUcpCN1Ov3HfVj70tR9yzbtP59yXH5tZOO2WFNYRVfnsBnD3+4ETE4pJRHJCXVLTE48oP1RnEGGa2k0KE+4+NmWbvi0is5ySQnriwYP1phtJU7sNzQ+b2UVAv5mtBH4b+F5yYYlIHignpCdOBvWmG0lTuyWF3wJ+lmiSu+uBMeCDSQUlIilq0KC5w0dUUkhRnAyyTgotSwphnME/u/vrgY8mH5KIpGrtJp7Zc5Cfv2IDF5xxHF+/dysLhgbYO17ib5UTUhMvh3oo4+qjliUFd58EymY2kkI8IpKB+Op0OIxo7u+L5kBSSSE98XoW9da1SFO7bQp7gQfN7FZgX7zR3X87kahEJFVxffZQmPsonhhP4xTSM5GTNoV2k8I3wo+IzEKlcqOSQmYhFU78GWRdfdRWUnD368KymSeHTT9294nkwhKRNMVXqfEsqQN9UXJQ9VF64tJaT1QfmdnriFZF2wIYcLyZXezudyYXmoikpVJ9pJJCZnqt+uhTwJvc/ccAZnYy8BXgjKQCE5H0xAOn4jaFAbUppC6uPprIe++jYDBOCADu/hOiNZdFZBY43NAcnRIG1Psodb1WUthoZn8H/EO4/25gYzIhiUjaSpUuqVFJoT9uU8j2/FQopV4ZvBa8H7iUaHoLgP8L/E0iEYlI6ibKKilkbbIclxR6oPdR2O8z7v5pqIxyHkosKhFJVXyVGjc0V9oUMouoeCbK+ag+ardNYQMwt+r+XOA73Q9HRLIwMWXwWlxSUENzevJSfdRuUhh2973xnXB7XjIhiUjaNHgte4cbmnuj99E+Mzs9vmNmq4ADyYQkImmbOs2FBq+l73CX1N5oaP4A8I9m9p/h/rHAu5IJSUTSFp+IKg3N/SoppC0v6ym0mxROBE4DTgDeDvw8aoMSmTXiaZvjLqlqU0jf4fUUeqP66GPuvgdYCLyeqDvqNYlFJSKpmlpSqLQpqKiQmjgxZz33UbtJYTL8/hXgb939n4E5yYQkImmLr07nVMYpxG0KmYVUOHEyyHqW1HaTwjYz+zxRO8ItZjbUwXNFJOfiE9Jgfx99Vt2moKyQll6b5uKdwJuBT7r7bjM7FlibXFgikqZS2bln6P3M//MxHh8CHoXPDBONULprKazdlHGEs1/c+yjr6qN211PYT9UiO+7+FPBUUkGJFM6VK2Hf9iO3z0/nhDwxWWbUxuo/WC8u6bpSTsYptFtSEJEkNTrxpnRCLmV8IiqEFol/orLyWm80NIvILDah6VCT1yLx52XlNSUFEVFJIQd6bZoLEZnFsr46lcMNzbO2+sjMjjez283sR2b2sJl9IGxfZGa3mtmm8PuYsN3M7LNmttnMHqiea0lEkjVRdnYyUv/B+UvTDaagilB9VAL+l7u/DDgTuNTMXgZcDmxw95VEHd4uD/ufC6wMP2vQiGkpkkYn3pROyKXJMm+Z80VYN1b5WXHwej591j3qjpqSeHxC2Q8vuJOFxHofVXdbdffnzewRYDmwGnhd2O064A7gw2H7lzyabOUuM1toZseG1xGZ3cKJ96M3PsiX7/4PHlz3Jo4aTm8Z9NKkMzhgNdv6THMfddX8pY17HxGNFemzKClMTJbp7+tPOcBIKm0KZraCaEK9u4FlVSf6p4Fl4fZy4Mmqp20N20QKY7Iy/026J+OJsjPYV3s66DPTiOZuadEd1d2ZLDtzw4SEWY5qTjwpmNkC4J+AD4ZJ9SpCqaCjb52ZrTGzjWa2cceOHV2MVCR7lUnRUq4+KE2WK1NbxKKkkGoYs1er7qjhDz13TlR5k2UPpESTgpkNEiWEL7t7PCL6mTBNBuF3/NfaBhxf9fTjwrYa7r7e3Ve5+6rR0dHkghfJQNzImHad8sSkVybBi5lp7qO0xCXDuXOiz2BWlhTMzIBrgUfc/dNVD90MXBxuXwzcVLX9vaEX0pnAmNoTpGhKGS3eXiqXGaxTUlBOSEc8eHDeYFxSyC4pJDnNxVnAe4AHzez+sO0jwMeBr5nZJcATRJPtAdwCnAdsBvYD70swNpFciksIR5QUEp4bqTTpDPRPbVPQegppOVxSiNsUZmfvo38DrMHD59TZ34FLk4pHpBc0bFNIeG6kiclyZbW1mNoU0hNXG+ahoVkT4onkSHxyKKU8F1GpqudLTG0KXdSoOyrAuhGWAluGgf8Ehjk8SiulWXKrKSmI5Egpoy6ppckyA8PhdBCqqh4AuC/8QCYnqFkj/N0uvf4+7nh0Ow/3vau952UwbbnmPhLJkYZtCgmr6X2U8TTes9meAxMsXjCUdRhNKSmI5MjhNoXsex9J9+05WGLJgnwvb6+kIJIjlTaFqdVHCc+NVK/3kXTf8z1QUlCbgkiONKw+CnXSn//uY/z5tx7lK79xJq/+mcVde9+JcpnBvi6VFDJeWjTP9hyc6LyksG7K7LUJ/x2VFERypDJ4rUGbQvx4t+fcj0oKbSSF+ATV7MTUrE1i3Uhhk4O7s+dAiZG5c9jhI43XxG4l4bYdlRdFcuRwSaH+ST/uvz4+MdnV952orj5qp0oqPsFfubLzNytog/V4qcyhyTJHzx3g7PJ6/uyV/86e/kV1931+oP72NCgpiOTIRKM2hSmPj5e6XFKorj5auwnWjfFXr93Y+okFPcFPx54DEwAcPTzIYL8xMekcPflc3X2PKtXfngZVH4nkSKsuqXGyONStpBDq/+8H+EH4AZi/lCW/uKE77yFA1J4AcPTcQeYM9Hc9sXeLSgoiOdKqTSGeE6drJ5Qm9f+jOe8l02vGDpQAOHp4gDn9lulUFs2opCCSI63aFCqLu5e626ZQz+hR00wKzaZ0KKJQGjuDMJXFV+B7wJ5Hj2n+vBYrtSVFSUEkRyYmm09z0ZWSQqMuo1Occu0Lp/f6azexbfcB5lz14vo9bFJadzo3Gvytj57c1fx5azdx5b88ytW3P8aas0/iI+e9NIHgjqSkIJIjcQmh0cprXWlo7uZVfIMT/DN7DvL28Wtqtt182Vm84riF3XvvWaBR19QdPsKCQ5Ps2h+1Q+wdL6UWk5KCSI60Wo4zHvHctYbmmWow3mD7nvEjto2F3jdy2CurEufRwwM8MP+yqD3HxuCKRVwBXDEMex46Bt62JZWYlBREcqTSptCgETJugB5PoU1hWkLV1JsJ9efVvkxhB641smX4osrtZ1kI+3bX3a9lVVMXKSmI5EjclpBZSWFdqMqYOrVCPfWqjlpVTcWD3uLnK0FULKZ+QkibuqSK5EipRZtCqdtdUqdhxcHr+Y/femrmJ/Si9FDqsYZ1lRREcqJc9sryl40Grx3qRkNzi1XAmL+04T6ORVUef1XnNXXVX1/4u6z60+/wSy9dyscffG3GATWnkoJITkxWLX3ZaGBTV0Y0r93EoeEljR/ftz06kdW5wjUaLP5TlKv+aSpNlnl23zhLj57a0NKB6c411SElBZGcqB6b0HCai3JcUphZQ/Ocgztb76QTfdfs3HsId1jaakBgq6qmFD4TVR+J5ET1amuNxylk36bQlEYz1wq9sV5A6I317ca77h1cxIK4Cq6dhv6EKCmI5ER16aDUqPqo3EGbQrPFbpIQ2iMmPraLlR/9Fh9648n89n3nNU8SV66c3W0RTY59xcHrOe2EhZx50mKuueMxfufsk7ksxdAaUfWRSE5Ulw4alhRKHZQUmi12k5R929m17xAAi+bPaa+LaoEtnDvIvMF+AI4aHsw4moiSgkhOtNOmMFFOYUTzDEsSz+2vSgrS1DHz5jB3TpwU8lFxk48oRKSmTWGiwYR4h8cpzKyhee/gIhZM1FnIpbpr6TTbB54LJYVj5s3ipNCFdajvGXo/o4+MwSPw68PAzeGHFsuiJrykqZKCSE4Wmq8uHTScOrtLI5o/+fJv8o37tvLAul9uvNPUY2+z8TNOCovbXaA+b+0K7XwfulA113iN5gbdfqf5Pp1S9ZFIFnXvdbTVplDuTu+j/YdKzB9K5ppwV6clhby1K8z0+1A9lqBBVdz+wcUdBpUelRREcqK6TaHlGs0TbVQfNVmkZd+hSeaFuuxue8+/nMJ7hoFPzfCFclKC61h1zGs34e685GPf5uLXrOAjD78V9m1n3sSz2cXXgkoKIjlR3abQco3mel1Wr1wZVfHEP01OqPvHu1tS+P7FP22+w7pGVSVNZFGCS2DE8HP7DjFeKnPsyHD+SkV1qKRQFL161VUgNeMUGrQpVC+y4+6YVTVKtnPCCft0u6Twn7sPdO21MtXqb9jJoLLwP7eYMHDt1pkElh6VFIoiJ/Xm0thEG9VHcVuDe+N2h3bsP1Ri3pwOrwmbdFU9/6aXTTuWWWva/1steh9BojOvJlZSMLMvAG8Btrv7z4Vti4AbgBXAFuCd7r7LosudzwDnAfuBX3P3+5KKTaRGRgukTzXZoqHZ3ZksOwuGBrid32DwT6ZRJRPsn05JYZq9kSq6+XfOW4+lqdZNY9nRJlVsZ16xgbNPXsJfXHDKDIJqT5LVR/8b+GvgS1XbLgc2uPvHzezycP/DwLnAyvDz88A14bdI8tZu4p4tz/GOz/07bz9tOZ9+16lH7pNC9VvLNoVPnsyW4e6U7PaPTzK/05LCTIW/00937uP1n7yDq951Cm877bj6+7aq22+0WE9uqkmnX4qrZ+G8QXbvT2c508S+Fe5+p5mtmLJ5NfC6cPs64A6ipLAa+JK7O3CXmS00s2Pd/amk4hOpFvf7f77RAukpVL/FiWCgz+pOnW3deq91I9wF8FD4gWRPmlNKAotCV9Xn9jU5yXVyrPG+jRJCp6+XhRalpaPnDqa2xnXaDc3Lqk70TwPLwu3lwJNV+20N245ICma2BlgDcMIJJyQXqRRKPEL4+YPZLS4ftyMMDfQ17H2UmCROmg2qQ44aHqC/zyrjGbqiWULoRFqzvHbYG2tk7iBPPrc/oWBqZdb7yN3dzDr+5rv7emA9wKpVq1L+z+lhOak3z6u4pLC3UUkhBXE7wtBg/4wakXOjQb1/36dO5rE52+HfiX5iMymtdOtEvnYTjz69h8XX/FyTEcczNI3/uZG5gzw0S0sKz8TVQmZ2LBB/ktuA46v2Oy5sk24J/2xX/sujXH37Y/z+r7yUX3/tSRkHlR/xCOG9B1NKCnWubN8I3DM0wuqBLzbskporra6qO63Kqa4GytCO58d5SVIJAaaV+BbO4uqjm4GLgY+H3zdVbb/MzL5K1MA8Vuj2hAQbyw5ORCebPSl9wXrFeLOSQhInqQYnxlEbi0oKDbqk5sraTd1fDCbDxWViO/eOJ/fi0yyZj8wdZP+hSQ6VyswZSHYkQZJdUr9C1Ki8xMy2An9IlAy+ZmaXAE8A7wy730LUHXUzUZfU9yUVV09odiU1w4RxIEyPsCetK+IeUWlorvd3aXY1nED1W6M2hdLcUQYO7Oj6+81Ip3XwWZUCOvicdjzfvaSww0e44Rc3cNkbZnbcI/OitRbGDkww2mpJzxlKsvfRf2vw0Dl19nXg0qRimVVm2LviYJwU8lhSyLA7YZwUxkvlzq7GEohreLC/blXB47/2A06+pkEXzpmabnJrVVqY2q6QVS+guAtrs+9S+P6t6eLbjtoYl925Cu5k+t/jK1fy3n3bee/U+aQS+r/QNBcFMx5XH1X3sslL3+5mCS/hOeSrZx3dO15i0UB2awEMDfTVbVOYmCyzw0e63wA6nXmJ2pW3rqCNxjfEjyX93t18XrNjmQFNc1EwleqjA1XVJL3StzvBeKrXJ0itsbmB4cF+Juu0KZQmnVeOX8NtF/6ksxN5q31nWqXTqz3Y4pNq/NPLuvi/oZJCXnSrn3ULleqjDPvjz1gCJZvqlcyeH0/hb9OgLn6HjzA82FdZN6FaXHoY6OvytdxMv3dJNDjXEyefvF2szDLFTgqtTsTtnmS6cZJK6Yt+IM9tCu1KoGTzmxvP5XeHw/KU66semL80mTEe4Xtx1sdvY9vuA9z6P8/mzk07+ZNv/oi3DvTXbWg+VAojnvvbmDBtaoxpDcpKUvy/1OtX9TlX7KTQ6p+k3X+iPFS/tHmCqnRJVe+jGgtKddYrhlDFMMa3H3qK//EP93Hey1/ALQ8+zaY/O5fB/plfscfJec/BEpPlcrRu749Ddc+6qh3nL6V0/vcADr9vs2Q19WIkrav5VqabnHq1iqoHFTspzAI/vOQJTjm+zRkZr1zJt/Zth+Fwf10bz6l3IkmqwTfHV7PxZGT/ZWQuEM0yOjJ3Zklhsuzcxm8wOjwGX4QzoPGsyfu2V8YuVJJCnmYJbbc0NZ2r/altInn9nrSKa7qJLeXjVVLocR21DXTri5XUFzScMJ5d90IWs7v+Phld7e4OV/QfeuAt/P7ws/CJKTvUS5QtqhX3HJjoqCfRa258NfBZBvo6qD5KS6cJaiYnurWbUmuDA2qTUhY99do53i6WpJQUely3hr7/0Rnf44v/bwu/euYJ/On5L482ZnQCfp2vr5mtdMvwRcm8UQcnll37DzGnv495hxqsrVvvdVpUK44dmOCYtt49MjS+E2D61VZ5mv9qpif2rNoXsiqdtar+62JcSgp50KRLYNmNvibzBtZ0LZ2BP7z3Nfzm0Ajr9v2frrzedJXLzt5D0zymTk9uHZyQxvZPsHDeIHSrfX7dCCum8bQtwxfB58KdTq9O81TdBLXxNLsCl1QVOyl0qw5wpldgTWJolhB2+EhXJ8katTF27e/idMaNND0BOD8dmsZUDvX64nejiiF8hrvjpJDgOK+O5bFefbqmm7Ca/e/Npr8PpFbSK3ZSCF/EN//lnTz69POYRWvffvqdp/D208N0Ah2cWHb4CH+w8kau+dUzWu/chRPWK8ev4f1dHm/wXKdz3Hdax9piIZSu9sCfyd93SpLZfeAQC+fOaZ4U0qznlkizZNLFi4JcSKmkV+ykEMSzInq4KN833sZo3zpGbezIWTYTPFFsGb4I7iL6qdbghFyeN0rf/uZX4TVL/jW7Mml1XHnoptsN4Ti/2s6+vXZss910u+EmOe1HDyh8UihNlnl2ytVxwyUZ27Bv6nOzOFHEs6lOSQzb//tDvOCqZQ2eFNm1/xDujpnB2k1cccsjrL/zcd78sy/gc++pKgG1889WvU/S8yjVxGPMdI3cyt9AJ/piyVPJICOFTwrP7jtUKSHEZjL3zd/tuAjWNehOmabqybKCF7TxtPFSmQMTk8wLi7rvDNMIz3hajDrxJGdmCWGHjzDnYImRuYNdiqcD8VVqHgaazQbttBvmrQE+Y8VNCqFaYBmwZfjw5h0+wl+P3zLtl13kOUgIM7Bl+CK44vD9P+g7hm9wdWqrPmVlxcHredtpy3ndi0f5wFfvZ8PecUb++mXpBxLPBtsOXdW2phN+x4qbFJqsfDWT6qPZZmF5F1uGL+LZ5xYCTyTbmDp/KQcmJplbZyzAxPASVu7+LF983yt5/YuXcvqf3Mp9kxd07a23DF8EjwCPwOph4OquvXTn2vn76gpXElLcpNBETbvAbOzaNg2VEcZJ/i32bWcuUWlt7kceZ8HQADvXncASxhg8uDM6cX8l2nVjw/kgCqDgDaGSLCWFOj7/2DntzQtUNCnVc4/aGPz5YgCWNNinb4btBrmXp9HHUijFTApZrRMr0i5VDUlGirnymqqDRETqKmZSEMmFBu0iqiKSDBWz+kgkD/IwnkVkimIlhV6Ym0a9nWYXNRhLjylWUsj7yba673mznj7rxvI74jXPsaVNYwmkBxUrKeRdJyeQ6ZQoqvu3Nztxz3Qd3TRKO6366ieVmFodm8YQSI9TUsiLqdUJraodpiaQVifBTqorql+73ZNr9ckwfn6zabVnkjTaOZZuzLPf6ASvBWFkFlNSyEqrK8puVjvUe6806rqbHcNMqsfa+dskOc++qoRkFlNSmKm43rjX6tHbPbFl1fCdZNLSMpAiDRUrKczkBNeq0XC29hqa6QLrjXRaPZYUXfWL1ChWUtv9AZwAAAa5SURBVJjOFX27PUg6WXoyiavQNK6su3k8OhmL5FKxkkIr3e45kuaJL4330olcZNbTNBciIlKRq6RgZm82sx+b2WYzuzyRN2lU1aGGRRGR/FQfmVk/0XpXbwS2AveY2c3u/qOuvpGqQEREGspTSeFVwGZ3f9zdDwFfBVZnHJOISKHkKSksB56sur81bKthZmvMbKOZbdyxY0dqwYmIFEGekkJb3H29u69y91Wjo6NZhyMiMqvkKSlsA46vun9c2CYiIinJU1K4B1hpZiea2RzgQuDmjGMSESkUc/esY6gws/OAvwT6gS+4+5+12H8H8MQ0324JsHOaz80bHUs+6VjySccCL3T3uvXvuUoKaTKzje6+Kus4ukHHkk86lnzSsTSXp+ojERHJmJKCiIhUFDkprM86gC7SseSTjiWfdCxNFLZNQUREjlTkkoKIiEyhpCAiIhWFTAqpTNGdIDPbYmYPmtn9ZrYxbFtkZrea2abw+5is46zHzL5gZtvN7KGqbXVjt8hnw+f0gJmdnl3kR2pwLOvMbFv4bO4PY2/ix34vHMuPzeyXs4n6SGZ2vJndbmY/MrOHzewDYXvPfS5NjqUXP5dhM/u+mf0wHMsfhe0nmtndIeYbwmBfzGwo3N8cHl8xrTd290L9EA2Meww4CZgD/BB4WdZxdXgMW4AlU7b9BXB5uH058Ims42wQ+9nA6cBDrWIHzgO+BRhwJnB31vG3cSzrgN+ps+/LwndtCDgxfAf7sz6GENuxwOnh9lHAT0K8Pfe5NDmWXvxcDFgQbg8Cd4e/99eAC8P2zwHvD7d/E/hcuH0hcMN03reIJYXZOkX3auC6cPs64PwMY2nI3e8EnpuyuVHsq4EveeQuYKGZHZtOpK01OJZGVgNfdfdxd/8psJnou5g5d3/K3e8Lt58HHiGaobjnPpcmx9JInj8Xd/e94e5g+HHgDcDXw/apn0v8eX0dOMfMrNP3LWJSaGuK7pxz4F/N7F4zWxO2LXP3p8Ltp4Fl2YQ2LY1i79XP6rJQrfKFqmq8njiWUOVwGtFVaU9/LlOOBXrwczGzfjO7H9gO3EpUktnt7qWwS3W8lWMJj48Bizt9zyImhdngF9z9dOBc4FIzO7v6QY/Kjz3Z17iXYw+uAX4GOBV4CvhUtuG0z8wWAP8EfNDd91Q/1mufS51j6cnPxd0n3f1UolmjXwW8JOn3LGJS6Pkput19W/i9HbiR6MvyTFyED7+3ZxdhxxrF3nOflbs/E/6Ry8DfcrgqItfHYmaDRCfRL7v7N8Lmnvxc6h1Lr34uMXffDdwOvJqoui5eSrk63sqxhMdHgGc7fa8iJoWenqLbzOab2VHxbeBNwENEx3Bx2O1i4KZsIpyWRrHfDLw39HY5Exirqs7IpSl1628j+mwgOpYLQw+RE4GVwPfTjq+eUO98LfCIu3+66qGe+1waHUuPfi6jZrYw3J5LtH79I0TJ4YKw29TPJf68LgBuCyW8zmTdwp7FD1HviZ8Q1c99NOt4Ooz9JKLeEj8EHo7jJ6o73ABsAr4DLMo61gbxf4Wo+D5BVB96SaPYiXpfXB0+pweBVVnH38ax/H2I9YHwT3ps1f4fDcfyY+DcrOOviusXiKqGHgDuDz/n9eLn0uRYevFzeQXwgxDzQ8AfhO0nESWuzcA/AkNh+3C4vzk8ftJ03lfTXIiISEURq49ERKQBJQUREalQUhARkQolBRERqVBSEBGRCiUFkYyY2evM7JtZxyFSTUlBREQqlBREWjCzXw3z2t9vZp8Pk5TtNbOrwjz3G8xsNOx7qpndFSZeu7FqDYIXmdl3wtz495nZz4SXX2BmXzezR83sy9OZ1VKkm5QURJows5cC7wLO8mhiskng3cB8YKO7/yzwXeAPw1O+BHzY3V9BNII23v5l4Gp3PwV4DdFIaIhm8fwg0bz+JwFnJX5QIk0MtN5FpNDOAc4A7gkX8XOJJoYrAzeEff4B+IaZjQAL3f27Yft1wD+GuaqWu/uNAO5+ECC83vfdfWu4fz+wAvi35A9LpD4lBZHmDLjO3X+vZqPZx6bsN935Ysarbk+i/0nJmKqPRJrbAFxgZkuhsm7xC4n+d+KZKi8C/s3dx4BdZvbasP09wHc9WgFsq5mdH15jyMzmpXoUIm3SVYlIE+7+IzP7faKV7vqIZkS9FNgHvCo8tp2o3QGiqYs/F076jwPvC9vfA3zezP44vMY7UjwMkbZpllSRaTCzve6+IOs4RLpN1UciIlKhkoKIiFSopCAiIhVKCiIiUqGkICIiFUoKIiJSoaQgIiIV/x9nXwPVbPhh1gAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A8YZpNIpmmtb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}